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NUMERICAL SOLUTIONS OF KNUDSEN FLOW 
ENTERING A CIRCULAR TUBE THROUGH 
A SMALL AXIAL ORIFICE 


SUMMARY 

The problem of determining Knudsen flow flux distributions in certain models 
pertinent to electric rocket thrustor design is discussed. A basic model of cir- 
cular cross section is established in which particles are assumed to enter one 
end through a small orifice and leave the other end, which is wholly open. 

Analyses are given for the cases of diffuse and specular reflections from 
the walls. For diffuse reflection, the resulting integral equations are solved 
by an iterative procedure on an IBM 7094 computer. Flux distributions are pre- 
sented for the rear wall, the side wall, and the open end of the basic model. 

For the case of specular reflection, approximate solutions are obtained by a 
method of ray tracing. The specular-reflection analysis is extended to complex 
models that closely approximate five variations of an electron-bombardment ion 
thrustor. Flux distributions for these models are found to be relatively insen- 
sitive to changes in propellant-distributor geometry. 


INTRODUCTION 

The subject of rarefied gas dynamics is receiving greater attention in re- 
cent years, largely because of the intense interest in exploration of outer 
space. In the study of electric propulsion, for example, a knowledge of the ar- 
rival rate per unit area, or flux distribution, of rarefied gases is of value in 
predicting thrustor performance. In general, a limited number of solutions to 
flux- distribution problems is available in the literature. Solutions are 
largely dependent on the assumed mode of surface reflection and on the geometry 
of the model under consideration. 

The study of free-molecule flow, or Knudsen flow, in tubes of right-circular 
cross section has perhaps received the most attention and dates back to the early 
work of Knudsen and Clausing. More recently several others have investigated 
this problem. Solutions to the "infinite tube" problem are available in most 
texts on kinetic theory (e.g. , ref. l). Solutions of efflux patterns from short 
tubes are presented in reference 2. In reference 3, perhaps the most comprehen- 
sive work of recent date on the Knudsen flow problem, the tube flow problem, 
among many others, is discussed. Other investigators have studied special varia- 
tions of the tube flow problem. For example, in reference 4, the problem of de- 
termining the pumping speed of a circular aperture in a diaphragm across an in- 




f ini te tube is studied. Recently, solutions for molecular-flow flux distribu- 
tions in right-circular tubes (ref. 5) and in tapered tubes (ref. 6) have been 
presented. 

Particle behavior in the free-molecule flow regime is of interest in the 
field of electric propulsion because, frequently, propellant-mass flow rates in 
electric rocket thrustors are such that particle densities are very low. Deter- 
mining the flux distributions of interest, however, is seldom a straightforward 
matter. An electrostatic ion thrustor that exemplifies the problem is shown in 
figure 1. Details of its operation and performance were first reported in refer- 
ence 7 and will be mentioned here briefly for orientation purposes only. Neutral 
propellant vapor passes from the vaporizer, through an orifice, and into the ion 
chamber. The propellant is ionized and a plasma is formed. By virtue of a po- 
tential difference maintained between the screen and the accelerator grid, ions 
are extracted, accelerated, and exhausted. 

The purpose of this report is to determine the neutral-particle flux distri- 
butions of certain axisymmetric models that arise in connection with this ion- 
thrustor configuration. Two analyses are made of a model consisting of a right- 
circular cylinder in which particles enter at one end through a small orifice on 
the axis and leave the downstream end, which is wholly open. The difference be- 
tween the two analyses stems from the mode of surface reflection assumed: dif- 

fuse or specular. Solutions of flux distributions are obtained on an IBM 7094 
computer. Because of its relative simplicity, the specular-reflection analysis 
is extended to models that are representative of the ion-thrustor configuration 
discussed. 


MODEL AND ASSUMPTIONS 

The basic model from which the two analyses are developed is shown in 
sketch (a). 



(a) 

Particles enter the right-circular, cylindrical enclosure through a small, 
thin. orifice located on the axis at the rear wall. The source of the entering 
particles is considered to be a large chamber in which the gas is in equilibrium. 
The entering particle flux is assumed uniform across the orifice opening. The 
pressure downstream of the open end is assumed zero except for that contribution 
that arises from the leaving particles. In one analysis, reflections from the 
walls are assumed diffuse; in the other, specular. In both analyses, the enter- 
ing flux from the orifice is assumed to follow the cosine law distribution and 
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the particle arrival rate at the vails is assumed equal to the leaving rate. The 
particle mean free path is assumed sufficiently large to ensure free-molecule 
flov everywhere vithin the model. 

Full consideration of the exact nature of the surface-reflection process is 
"beyond the scope of this report. Reference 8 suggests that the occurrence of 
specular reflection (wherein the angle of reflection of a particle leaving a sur- 
face equals the angle of incidence) is unlikely. Reference 9, however, points 
out that specular reflection may he of practical importance in very high speed 
flows. It is generally agreed among experts that diffuse reflection (wherein the 
direction of reflection of a particle from the surface is unrelated to the direc- 
tion of incidence) is the more probable mode. 


ANALYSES AND PROCEDURES 


Fundamental Relations 

The analyses to be discussed are based on fundamental relations developed 
from elementary kinetic theory. 

The rate of particle flow from a differential area dSj_ in a direction 0-j_ 
arriving at an area dSj is 

_ _ 

d( n i j) = n-^ cos Si dS ± — (1) 

where Vj_ is the isotropic mean particle speed, n^ the particle density, and 




cos dSj 



(2) 


(All symbols are defined in appendix A. ) The subscript notation is described in 
sketch (b). 
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The arrival rate, or flux, at dSj from dS 4 is, from equations (l) and 


With n^ 


d(n., n -) 

’Ssf- 


(l/4)n i v^, equation (3) 


= i cos cos Sj 

can also "be written as 


dS ± 



. • . n-; dS-? 

d(n j) = ~ cos ®i cos 0 j ~ 2 ~ 

L ij 


( 3 ) 


(4) 


From equation (4), the flux through a differential area dS^ of a hemisphere of 
radius 2-^, from a small orifice centrally located on the "base, can he written as 


fi4 Ai n-t 

d(n h ) = — ri cos Q ± = — - cos Q ± (5) 

jti h ni h 


In equation (5), the cos = 1, and the orifice area A 4 is considered small 

2 

enough to he almost a point source, that is, Aj_ « kI-^. 

Equation (4) is the basic relation from which the subsequent analyses are 
developed. It will he used as written for the diffuse-reflection analysis and as 
rewritten in the form of equation (5) for the specular-reflection analysis. 


Diffuse -Reflection Analysis 

In formulating the problem described herein, a set of two simultaneous in- 
tegral equations arises from considering the contributions from all sources to 
the local arrival rate (assumed equal to the leaving rate) on each surface. 

The equations that describe these local fluxes are given next in general 
form and are derived in detail in appendix B. Model surfaces and pertinent var- 
iables are identified in figure 2. The local flux on the rear wall n 4 (r 4 ) in- 
cludes contributions from all local points on the side wall and can be expressed 
as 

n 4 (r 4 ) = f n 3 (x)F(x,r 4 )dx (6) 

*0 

(rear (side wall) 

wall) 

On the side wall the local flux n 3 (x 3 ) includes contributions from the orifice, 
the rear wall, and the side wall Itself: 
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n 3 (x)I(x, x 3 )dx 
(side vail) 


(7) 


h 3 (x 3 ) 


( side 
wall) 


hiG(x 3 ) + 

(ori- 
fice ) 


f n 4 (r)H(r,x 3 )±r + f 
\/2 J 0 
(rear wall) 


The local flux in the plane of the open end hg(sg) includes contributions from 
the orifice, the rear wall, and the side wall: 


n 2 (s 2 ) 

( open 
end) 


niJ(s 2 ) + 

(ori- 

fice) 


f n 4 (r)K(r, s 2 )dr + f n 3 (x)M(x, s 2 )dx 

*4-l/ 2 *t) 

(rear wall) (side wall) 


(8) 


In these equations, the functions F, G, H, I, J, K, and M are the kernel 
functions of the model geometry (see appendix B). The side wall is the wall of 
the cylinder. The flux entering through the orifice, rrj_, is assumed known, that 
is, a constant. 


Equations (6) and (7) must he satisfied in order to obtain a solution to 
equation (8). The basic difficulty encountered in seeking solutions to equa- 
tions (6) and (l) stems from the fact that the functions h 3 (x) and h 4 (r) - the 
flux distributions along the side and rear walls, respectively - are unknown. 
Thus, direct integration is not possible. 


Described in reference 10 is a numerical method of solution to the problem 
of determining the thermal-radiation characteristics of cylindrical cavities. 

The simultaneous integral equations for the apparent emissivity distribution that 
are given in reference 10 are similar in structure to equations (6) and (7). The 
procedure followed herein to obtain solutions to equations (6) and (7) employs 
techniques similar to those used in reference 10. Solutions are obtained by 
iterative numerical integration performed on the IBM 7094 computer of the NASA 
Lewis Research Center. 


Diffuse-Reflection Procedure 

For actual computation, equations (6) and (7) are first normalized by divid- 
ing through by n-j_, the flux through the orifice. To start the iterative method 

of solution, an initial guess is made of the flux distribution h 3 (x) in equa- 
tion (6). The significance of the value chosen for the initial guess of n 3 (x) 
will be discussed later in this section. The function n 4 (r 4 ) is then calculated 
by supplying values of r 4 , d^/2 < r 4 < R to the integrand function, for in- 
stance f(x, r 4 ) = n 3 (x)F(x, r 4 ), and then by numerically integrating f(x, r 4 ) by 
Simpson's rule to obtain pointwise values. 

To obtain values of n 3 (x 3 ) from equation (7), a similar procedure is fol- 
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lowed using the values of n^/r) produced from equation (6) and the initial guess 
of hg(x). Values of Xg, 0 < Xg < L are supplied in G, H, and X of equation 

(7). The function G is prescribed: the integrand of the first integral "be- 
comes a function, for instance, h(r, X3); and the integrand of the second integral 
"becomes a function, for example, i(x, x 3 ). The functions h(r, Xg) and i(x, Xg) 
are numerically integrated, added to G(xg), and values of hg(xg) are thus ob- 
tained. The newly tabulated rig(x) function is then used as the "next guess" in 

equation (6) and the cycle is repeated until convergence is obtained, that is, 
until the values of hg(x) obtained from equation (7) agree with the values of 

hg(x) supplied previously in equation (6). 

With the flux distributions, hg(x) along the side wall and h^(r) along the 
rear wall, determined, the values of the local flux distribution, hg( Sg), across 

the open end can then be determined by numerical integration of equation (8). No 
iteration is required in this case; however, the integration must be done once 
for each value of the radius, sg. 

For the iterative integration of equations (6) and (7), the r and x in- 
crements used were as follows: on the rear wall, l/42 of the cylinder radius; 

on the side wall, l/40 of a unit length; and, on the open end, l/40 of the cylin- 
der radius. The increment used on the rear wall was chosen because the Simpson's 
rule integration requires an even number of increments. In the first integral of 
equation (7), the r integration is carried out from d-^/2 to R. With 

d-j_/2R = l/21, 40 increments were available. The second integral of equation (7) 
was broken into either two or three parts, depending on whether the Xg point 
being calculated was numbered odd or even. When Xg was an odd-numbered point, 
an even number of increments was available to the left and right of Xg and 
Simpson's rule integration was applicable. When Xg was an even-numbered point, 

Simpson's rule integration was applied over the range from x = 0 to the x 
step value preceding the Xg point, and from the x value after the point out 

to x = L. Over the interval between, trapezoidal integration was applied. 

Because of the complex structure of equations (6) and (7), trial solutions 
were investigated in an attempt to establish convergence criteria and overall ac- 
curacy of the results. 

To determine the effect of the initial guess of the hg(x) function that 

must be supplied in equation (6) to start the iterative procedure, several dif- 
ferent functions were investigated. Functions that had positive slopes, negative 
slopes, and even a value of zero, were tried. In all cases, convergence to the 
same final answers was indicated. No stability problems were encountered; how- 
ever, the rate of convergence depended on the particular initial guess selected. 
The number of iterations required for the final solutions was reduced to low 
values (as indicated in table I) by noting that the constant- slope portion of the 
side-wall relative -distribution curves could be approximated by an initial guess 
of rig(x)/n-]_ « 0.00l{ [(L/R) + 1] - (x/R)). 
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The "ratio test" indicated in table I constitutes a check on the numerical 
results from the law of mass conservation. This check ratio is derived from a 
comparison of the flow in through the orifice and the flow out through the open 
end plus the small fraction returning hack through the orifice. Because of the 
area ratios, this latter fraction would he expected to he small. Expressed 
mathematically, the comparison is 


2it 



ii 4 (r)r 




( 9 ) 


or, in terms of relative flux distributions, 




( 10 ) 


Eor exact agreement, a ratio of the left side to the right side of equation (10) 
should equal unity. To make the desired comparison it is necessary to calculate 
the "return-flux" distribution across the orifice opening from equation (6) and 
then perform the integration indicated in equation (10). A trial calculation for 
a length to radius ratio (L/R) of 5 gave a return- flux ratio of 0.0046 and an 
open end ratio of 0.9713, which is a total of 0.9759. Although numerical inte- 
gration that involved several integrals was required to obtain this ratio, the 
result seems very good. Other trial solutions showed that this degree of agree- 
ment was not possible with coarser increments, while the use of finer increments 
resulted in little or no improvement and added greatly to machine computation 
time. Since the return-flux ratio was negligible by comparison and calculation 
also required additional computer time, it was neglected in obtaining the final 
solutions. 

As a result of the preliminary evaluations discussed, a convergence crite- 
rion was established for the final runs and consisted of the following two 
checks: 

(1) After 20 iterations of equations (6) and (7), equation (8) was solved. 
The flux distribution was then integrated as in equation (10), and a ratio of the 
left side to the right side was formed. If the ratio obtained was greater than 
an arbitrarily selected value - either 0. 96 or 0. 97 was used - then the iterative 
procedure was stopped and the results were listed. It was necessary to decrease 
the value of the comparison ratio as the length to radius ratio increased because 
of inherent inaccuracies of the numerical integration. 

(2) If the ratio obtained in item (l) was less than the preassigned value, 
an additional 20 iterations were performed on equations (6) and (7) and a new 
ratio comparison was made. If the condition of item (l) was satisfied, the re- 
sults were listed. If not, the ratio obtained was compared with the ratio of the 
previous 20 iterations. If the change in the ratios compared was not greater 
than 0.0001, iteration was stopped, and the results were listed. Otherwise an- 
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other 20 iterations were accomplished and new comparisons were made. 


As a precaution to prevent endless iteration if the ratio test failed, the 
program was to he stopped after 100 iterations and the results listed. In the 
final solutions that employed the initial guess given previously, this precaution 
was not needed. 


Specular-Reflection Analysis 

An approximate method of solution of the specular-reflection problem was 
developed that is readily adaptable to numerical techniques. 

If equation (5) is multiplied through by dS^, where, from sketch (c), 
dSfa = sin 0^ d0^ dcp, the result is 

_ n-t 

d( n h) = — cos 0^ sin 0^_ d0^ dcp (ll) 



In equation (ll), d(n^ 1 ) is the arrival rate (particles per unit time) on 

the differential element of the hemisphere. The arrival rate passing through a 
differential band is 


E -1 

d (“b) = — 

r 2rt 

/ cos 0j_ sin 0j_ d0-j_ dcp 


*u 


= *i 

2 cos 0j_ sin 0^ d0^_ 


= n ± 

sin 2 0j_ d©i 

(12) 
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For the analysis, it is assumed that the particles initially leaving the 
orifice in a direction near and about 6 ^ are restricted always to follow the 

direction traced by the "ray" associated with 0^. The arrival rate (particles 

per unit time) on an arbitrarily narrow band on the hemisphere (sketch (c)) is 
then written from equation (12) as 

AiTb ~ n^ sin 20^ A 0^ (13) 

To find the flux (particles per unit area per unit time) through an annulus, 
such as shown in sketch (d), it is necessary to sum the arrival-rate contribu- 

Side wall— \ 



tions from all rays passing through the annulus and to divide by the area of the 
annulus, that is, 


E^ 


(14) 


Substituting equation (13) into equation (14) and noting that 
Ap = 7t[(r + Ar) 2 - r 2 ] yield 




5=3 


The summation i through j 
An alternate approach to the 


n^ A 0 


± 


sin 20^ 


(15) 


rt Ar(2r + Ar) 
includes all rays passing through the annulus, 
specular- reflection problem was suggested by 
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Eli Reshotko of Lewis and is given in appendix C. Although the method yields 
an analytic expression for the flux distribution across the open end of the basic 
model, it is not readily extended to handle more complex models. Agreement of 
results between the two methods for the basic model, however, indicates that the 
method used herein can reasonably be extended to more complex geometries. 


Specular-Reflection Procedure 

The flux calculation which is largely a matter of bookkeeping, was accom- 
plished using the IBM 7094 computer. For the model shown in figure 2, a mesh was 
overlayed. The mesh size used in the r-direction prescribed the values for the 
deno mi nator of equation (15) for each annulus. The value of A 6 is prescribed 
by dividing the rays that emerge from the orifice and range in value from zero to 
at / 2 into an arbitrary number of rays such that the first ray is A9/2, the next 
is [(A0/2) + A9], the next is [(A0/2) + 2 A0], and so on. 

The procedure is begun by following the first ray throughout its entire path 
and noting the r-count location at each x-count and storing this information. 

The next ray is then followed and its path recorded in a similar manner. When a 
ray arrives at the side wall it is reflected, the angle of reflection being set 
equal to the angle of incidence. In the same manner, the rays, which are actu- 
ally representative of cones, are reflected from the axis. Thus, it is possible 
to accomplish the entire calculation by following rays in the r,x plane. For ex- 
ample, the r-count position of the initial ray at any x-count position (see 
sketch (e)) is given by r = x tan Q ■$_. 
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These position values are recorded as the x-count is increased in steps un- 
til the value r = R is reached. The equation of the first reflection is 
r = 2R - x tan 0j_, and the procedure of storing r-count locations of Q±, as the 
x-count is increased, is continued until the value of r = 0 is reached. The 
ray is reflected from the axis, the equation of the second reflection being 
r = x tan 0j_ - 2R. The equation of the third reflection is r = 4R - x tan 9±, 
and so on. With each reflection from the vail r = R, the sign of the x tan 9± 
component alternates from positive to negative, and the R component increases 
by 2R. With each reflection from the axis, the sign alone alternates. In this 
manner, the rays passing through a given r-count location are recorded, 

and after completing the tracing of all rays, equation (15) may be applied at the 
x-planes of interest. 

Flux distributions on the side vail can be obtained in a similar manner by 
counting and storing ray information in x-count locations at r = R and then 
dividing by the incremental cylinder surface area 2jtR Ax. 

Because of questions previously discussed concerning the plausibility of the 
specular reflection process as veil as recognition of the approximation of the 
method, the discussion on the accuracy of the method is confined to the details 
of the calculation itself. The assumption that rays initiate on the axis at the 
plane of the rear vail obviously excludes application of the method to problems 
other than those in vhich the orifice diameter is very small compared vith the 
cylinder diameter. 

The accuracy of the calculation is strongly dependent on the size of the A 8 
Increment and somevhat less dependent on the size of the Ar increment. The in- 
tegral over 0 in equation (12) is equal to unity and the summation over 0 
given in equation (15) approximates this Integral. For example, if R vere di- 
vided into 10 increments, vhile at the same time the A 0 increment vas selected 
to be 45° (tvo rays), there vould be several locations through vhich no ray vould 
pass. Also, if only tvo rays vere used (A0 = 45°, 0-j_ = 22.5°, 0g = 67.5°), the 

summation over 0 vould equal 1.112. Therefore, to include all locations and 
to obtain a good approximation of the Integral, it is desirable to choose A0 
as small as is practical. Optimization of the A 0 increment . and Ar increment 
vas investigated. It vas found that a Ar increment of from 20 to 25 units of 
R, used in conjunction vith a A 0 increment of 0.125° (i.e., 720 rays), gave 
good results - measured in terms of the degree of "smoothness" of the distribu- 
tion curves. With this choice of A 0, the summation over 0 equals unity to 
five significant figures. A A 0 increment of 0.0125° (i.e., 7200 rays) vas in- 
vestigated and vas found to have little effect on the results. 

It is apparent from inspection of both the diffuse- and specular- reflection 
equations that flux distributions are dependent only on the cylinder length to 
radius ratio. All results to be presented have thus been normalized by the cyl- 
inder radius, R. 
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RESULTS AND DISCUSSION 


Diffuse Reflection 

V al ues of the relative flux distribution along the rear -wall, the side wall, 
an d the open end of the model shown in figure 2 were obtained for cylinder length 
to radius ratios of 1, 2, 3, 4, 5, and 10. The orifice to cylinder radius ratio 
in all cases was l/21. Results are plotted in figures 3 to 5, and some tabulated 
values are given in tables I and II. These values are given to three significant 
figures, although machine computation gave eight. 

Turning first to the results shown in the figures, it is seen from figure 3 
that the distributions on the rear wall for various L/R are similar in form and 
differ pri mar ily in magnitude. The flux- distribution curves across the open end 
(fig. 4) approach the shape of the infinite-tube-distribution curve (ref. 1) as 
the L/R increases. In making this comparison, the infinite tube distribution 
from reference 1 has been normalized to the L/R = 10 value on the axis. The 
distributions in figure 4 are quite similar for length to radius ratios greater 
than two. Inspection of the distributions of figure 4 indicates that, as would 
be expected, the total relative arrival rates at the open end are nearly equal 
to unity for all L/R solutions (see table I). 

The flux-distribution curves along the side wall, given in figure 5, are all 
similar, each exhibiting an initial increase followed by a uniform decrease as 
x/R increases. As mentioned in the Diffuse-Reflection Procedure section, it was 
found that the constant- slope portion of these curves beyond the peak is closely 
approximated by the relation fi3(x)/n-|_ = 0.001{[(L/R) + 1] - (x/R)}. 

Because of the steep slope of the curves at the intersections of surfaces, 
end-point values of the curves of figures 3 to 5 are not shown in the figures but 
are given in table II. Values at the Intersections of surfaces, for example, 
fr^(R) and £13(0) are different. The mathematical reason for these differences 

cannot be rigorously explained from equations (6) and (7) but is somewhat appar- 
ent If the equations are compared, for example, for values of r^ = R and 
X3 = 0, respectively. At this point, the equations become similar in form, but 
each contains a different kernel function. In reference 10, it was found that 
similar differences in corner values occurred in solutions to the problem of ap- 
parent emissivity of a cylindrical cavity. The opinion expressed in reference 10 
is that the finding seems physically reasonable, based on the idea that an ob- 
server approaching the corner along the side wall would have a different view 
from that of an observer approaching the corner along the rear wall. A similar 
idea may apply to the problem herein. 

Extending the diffuse-reflection analysis to models such as that represented 
by the electron-bombardment thrustor shown in figure 1, wherein a distributor and 
screen grid are included, poses a formidable problem. The governing integral 
equations may be readily written; however, instead of two equations, three or 
more arise (some of which contain integrals having variable limits of integra- 
tion) that must be simultaneously satisfied. While solution by the iterative 
method is possible in principle, it may be much more difficult to accomplish 
practically because of the additional number of unknown distribution functions 
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involved. The iterative procedure can no longer he initiated with a guess of 
only one unknown distribution function. 


Specular Reflection 

Values of the relative flux distribution across the open end of the basic 
model of figure 2 are given in figure 6(a). These results may be compared with 
the results for diffuse reflection given in figure 4. Both sets of curves show 
maximums in the center (axis) and minimums near the wall; however, the portion 
of the curves between have opposite curvatures. The shape of the "specular" 
curves was nearly the same for all cases investigated. In this model, particles 
proceed in the forward direction only, so that the flux distribution on the rear 
wall is zero. 

Values of the open-end distribution determined from the approximate method 
(eq. (15)) are compared with those of the analytic method (eq. (C5)) in figure 
6(b) for an L/R of 5. In determining values from equation (C5), 100 terms were 
used (i.e., a = 100). It can be seen from figure 6(b) that the results from the 
approximate method compare quite well with the analytic results. Values for an 
L/R of 1 determined from equation (C5) were practically Identical with those 
determined from equation (15). As mentioned earlier, the analytic method is use- 
ful for comparative purposes but cannot be easily extended to more complex 
models. 

A side-wall distribution for an L/R of 4 is shown in figure 7 and compared 
with the diffuse-reflection distribution. The variation of the specular- 
reflection curve with x/R is typical of those obtained with other L/R values. 

This method of analysis was also applied to models representative of the 
electron-bombardment thrustor shown in figure 1. Ion-chamber pressures in this 
thrustor have been estimated to be of the order of 10”^ millimeter of mercury. 

The thrustor uses mercury vapor as a propellant. The mean free path of the mer- 
cury atom at this level of pressure with thermal equilibrium assumed is approxi- 
mately 25 centimeters compared with chamber dimensions of about 10 centimeters. 

It is likely, therefore, that molecular flow occurs In the chamber. Thermal en- 
vironment and the ionization processes, however, could have an appreciable effect 
on the flow. 

Various propellant "distributors" have been tested in this thrustor in an 
effort to alter neutral propellant distribution and improve thrustor performance. 
Five models that approximate several of these configurations, along with a sketch 
of the mesh structure used for calculations, are shown in figure 8. The radial 
flux distribution was of primary interest in these studies so that a coarse x- 
count mesh was used. In addition to the distributor, a "screen-grid" was located 
across the downstream end of the models. To preserve axial symmetry, the screen- 
grid was approximated as a series of concentric annuli. The reflection and re- 
reflection of rays was accomplished as previously discussed. Re-reflected rays 
are those that strike a downstream surface and return in a negative x-direction. 
Because of these rays, two modifications to the computer program were required to 
prevent the computer from continuously looping. In the first modification, rays 
that returned to the "orifice" by re-reflection were deleted. In the other modi- 
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fication, rays that became "trapped" by re-reflection were deleted after 200 re- 
reflections in the region between the rear wall and the distributor and after 
10 re-reflections in the region between the distributor and the screen. These 
limiting val ues were selected after considerable investigation of trial solutions 
run to much higher limits. Both of these "deletion" procedures may be justified 
to the extent that, under steady-state conditions, the particles, represented by 
the rays, must eventually find their way out of the chamber, otherwise, chamber 
pressures would become infinite. The number of deletions resulting from these 
modifications was actually a small percent of the total number of 720 rays fol- 
lowed. In the worst case, model E, 40 rays (5.6 percent of the total) were de- 
leted. Of these, only four were deleted in the region between the distributor 
and the screen. This region is of particular interest because the ionization 
process in the real thrustor takes place therein. The weighted value of the de- 
leted rays, however, must be considered to evaluate the effect of these deletions 
properly. A summation across the screen openings shows that, in the worst case 
(model E), a total weighted value of less than 4 percent was being deleted. 

Elux distributions at four axial stations for the five models studied are 
plotted in figure 9. Shown in figure 10 are the average flux distributions (two 
values per opening) through the screen apertures. The curves of figure 9 have 
the same general form as those of figure 6(a), indicating that the flux distri- 
butions are relatively insensitive to distributor geometry. In this respect, the 
specular- ref lection results are in agreement with thrustor test results reported 
in reference 11. Therein, only slight differences in performance were noted in 
comparing test results from configurations similar to models D and E. As can be 
seen from figure 8, these two models represent opposite concepts of propellant 
introduction. 


CONCLUDING BEMARKS 

The purpose of this report was to determine particle flux patterns that 
arise in connection with Knudsen flow in certain geometric configurations of in- 
terest in the field of electrostatic propulsion. A basic model of circular cross 
section was established in which particles were assumed to enter one end through 
a sttirI 1 orifice and leave the other end, which was wholly open. Two modes of 
surface reflection were studied: diffuse and specular. Numerical means were 

used to obtain solutions. The specular-reflection analysis was extended to more 
complex models that included a distributor and a screen. These models closely 
approximated variations of a Lewis Research Center electron-bombardment ion 
thrustor. Althrough the diffuse-reflection analysis could in principle be ex- 
tended to these models, the additional integral equations that arise make prac- 
tical solution of the problem extremely difficult. 

It was found that, with diffuse reflection, the shape of the flux- 
distribution curve across the open end of the model approached that of the dis- 
tribution curve for an infinite tube as the tube length to radius ratio in- 
creased. Eor the model studied, the relative flux distribution along the side 
wall was found to be linear in part and fitted a rule-of-thumb relation. Differ- 
ences in flux values at surface intersections were found and were postulated to 
be due to differing "view factors. " 
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Flux- distribution curves for both the specular- and diffuse-reflection anal- 
yses showed maxi mum values on the axis and minimum values at the wall: between, 

the distribution curves had opposite curvatures. Increasing the length of the 
model as well as inclusion of distributors with various openings had little ef- 
fect on the specular-reflection flux distributions. This insensitivity to dis- 
tributor geometry has also been noted in measured performance of electron- 
bombardment thrustors. Since diffuse reflection is generally agreed to be more 
likely than specular reflection, it must follow that flux distributions obtained 
from the diffuse-reflection analysis should be closer to the actual case. 


Lewis Research Center 

Rational Aeronautics and Space Administration 
Cleveland, Ohio, September 24, 1963 



APPENDIX A 


SYMBOLS 

A area, surface 

d diameter of orifice 

F function in eq. (6), defined in appendix B 

f distance, variable, used in derivations in appendix B 

G function in eq. (7), defined in appendix B 

g distance, variable, used in derivations in appendix B 

H, I functions in eq. (7), defined in appendix B 

J, K functions in eq. (S), defined in appendix B 

L length of model 

Z length, variable 

M function in eq. (8), defined in appendix B 

m distance, variable, used In derivations in appendix B 

N normal to surface 

n density, particles per unit volume 

n flow rate, particles per unit time 

n flux, particles per unit time per unit area 

P radius of model 

r radial distance, variable 

S area, differential 

s radial distance, variable, plane of open end of model 

v mean particle speed, distance per unit time 

x distance, variable 

a number, 1, 2, 3, . . . , «> 

(3 angle, variable, used in derivations in appendix B 
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6 angle, variable, between IT and l, measured from IT 

cp angle, variable 

cd solid angle 

Subscripts: 

b band on hemisphere 

h hemisphere 

i, j number 

p plane, normal to x-axis 

1 orifice 

2 open end 

3 side wall 


4 


rear wall 



APPENDIX B 


DERIVATION OE INTEGRAL EQUATIONS USED IN 
DIFFUSE- REFLECTION ANALYSIS 

Equations (6) to (8) are established by considering the contributions to 
the local flux on a surface from all possible sources and then by integrating 
over the contributing surfaces (i.e., by application of eq. (4)). Subscript num- 
bering in the following derivations and sketches is consistent with that on the 
model shown in figure 2. 

I. Equation (6) derivation. 

The contribution to the local flux on the rear wall at a point r^ from a 
local point on the side wall (see sketch (f)) is 


^ 4 ^ 4 ) = T C0S 6 43 C0S 0 34 ,2 


dS 3 


'43 


(Bl) 



dSj = R dp dx. 


cos 9 


x 


43 =T 


43 


cos 9 


R — r a cos p 


34 “ 


'43 


f2 _ 7 2 _ x 2 

x - - x 


Using these relations in equation (Bl) yields 
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r 4 cos dp ^ 


43 


'43 


Also, from the law of cosines, 

,2 


COS P ss 


R + r 4 - f c 

2 W A 


R 2 + r| + x 2 
2Rr, 


- I 


43 


(B2) 


from which 

l| 3 = R 2 + r| + x 2 - 2Rr 4 cos p 

The integral over all of the side wall yields the total local flux at a point on 
the rear wall: 


n . (r . ) = — 
4 4 jt 



R - r cos p 

2 dp dx (B3) 

(r 2 + r 2 + x 2 - 2Rr^ cos p) 


The flux distribution on the side wall is independent of p, and the p integra- 
tion can he performed (ref. 12) to obtain equation (6): 


n4(r 4 ) 



n 3 (x)F(x,r 4 )dx 


where 


F(x,r 4 ) 



A condition on the 


P integration given is that 


( 6 ) 


(B4) 





(-2Rr 4 ) 2 


This condition is everywhere satisfied except at the point x = 0, r 4 = R. 

II. Equation (7) derivation. 

The contributions to the local flux on the side wall at a point x 3 come 
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from three sources: the orifice, the rear wall, and the side wall itself. The 
equation -will he developed by considering each of these in order, with notation 
as shown in the sketches. 

A. The contribution from the orifice, considered as a point source, 
(see sketch (g)) is 


ni dSj_ 

ldh 5 (x 3 ) = — cos 0 31 cos 0 13 -g- 

l 15 


(B5) 



The following relations are established from sketch (g): 


/ 


ds i ■ i 4 


cos 0 


x. 


13 "T 


13 


cos 6 


R 


31 ~ 7 


13 


T 2 _ J>2 yZ 

4 13 ~ 3 


Using these relations in equation (B5) yields the total contribution from 
the orifice 


1 n 3 (x 3 ) = n 1 G(x 3 ) 


(B6) 


where 


G(x 3 ) - 


x 3 e 


4=(r 2 + x! 


ir 


(B7) 


B. The contribution from a local point on the rear wall ( see sketch 
(h)) is 
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4 <%(x 3 ) 


n 4 d-S 4 

= — cos e 34 cos e 43 — 



The following relations are established from sketch (h): 
dS 4 = r dp dr 


x 3 

cos e 43 = - 


cos 9 


34 


E - r cos P 
x 43 


43 


y 2 7 2 

D 6 43 3 


Using these relations in equation (B8) yields 




dp dr 
2 

43 


Also, from the law of cosines, 


cos p = 


r 2 4 p{2 _ g2 r + + x 3 - l 43 


2rR 


2rR 


(B8) 


(B9) 


from which 


+ R^ + Xg - 2rR cos p 


The integral over the rear wall yields that contribution to the total flux 
at a local point x 3 on the side wall: 


4 n 3 (x 3 ) 



R - r cos 3 


(r 2 + R 2 + x§ - 2rR cos p)* 


dp dr 


(BIO) 
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The 


0 integration can "be performed via reference 12 to yield 


/•* 

4 n 3 (x 3 ) = I n 4 (r)H(r,x 3 )dr 

\/Z 


where 


H(r, x 3 ) 



A condition on the 


0 integration in this case is that 


(r 2 + R 2 + x 2 ) 2 > (-2Rr) 2 


The condition is satisfied everywhere except at the point x 3 = 0, r ; 

C. The contribution from a local point (identified as point 5 in 
(i)) on the side wall itself (see sketch (i)) is 

n 3 dS 3 

5 dh 3 (x 3 ) = — cos 0 35 cos e 53 -g- 

z 53 



The following relations are established from sketch (i): 
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(Bll) 


(B12) 


: R. 

sketch 

(B13) 


Mill I I 



dSg = R dp dx 


cos G 35 = cos S 53 


R + m 
Z 53 


,2 « 2 / \ c 

h = 2 53 " (X 3 “ X) 


m = -R cos p 


Using these relations in equation (B13), and noting that n 5 (x) = n 3 (x)^ 
yields 

,. / X n 3 R 2 (l - cos p ) 2 R dp dx ,„,. N 

5<^3(x 3 ) = — I~ 2 if ( B14 ) 


'53 


'53 


from the law of cosines 


cos p = 


R 2 + R 2 - h 2 
2RR 


from which 


h 2 = 2R 2 ( 1 - cos p) = - (xj - x ) 2 


or 


l | 3 = 2R 2 ( 1 - cos p) + (x 3 - x ) 2 

The integral over the side wall yields that contribution to the total local 
flux at a point x 3 on the side wall: 


5^3( x 3) = “ 


/ L 

n 3 (4)R 3 I 

J n 


(l - COS P ) 2 

^(x^ - x ) 2 + 2R 2 (l - cos 


«T 


dp dx (B15) 


After some manipulation, the p integration can be performed via refer- 
ence 12 to yield 


where 


5 h 3( x 3) 



n 3 (x)l(x, x 3 )dx 


I(x,x 3 ) 



I x 3 - x l 


(x 3 - x ) 2 + 6 R 2 

£(x 3 - x) 2 + 4R^J ^ J 


(B16) 


(B17) 
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1 1 II ■■■■■■■Ill ■ 


and where the magnitude sign preserves the positive sense associated with 
the leg of the triangle in sketch (i) and permits the contribution over the 
entire surface of the side wall to he written under one integral. A condi- 
tion on the (3 integration in this case is that 


J(x 3 - x) 2 + 2R 2 J 2 > (-2R 2 ) 2 

This condition is satisfied everywhere except at the point x 3 = x. At this 
point, the latter half of the term in brackets in equation (B17) becomes 
zero. Rote, however, that this portion of the bracketed term is always pos- 
itive and approaches zero as x approaches x 3 . Finally, combining equa- 
tions (B6), (Bll), and (B16) yields equation (7): 

/* 

h 3 (x 3 ) = fi^Xg) + / n 4 

*V 2 


(r)H(r, x 3 )dr 




n 3 (x)l(x, x-^ ) dx 


(7) 


III. Equation (8) derivation. 


The derivation of equation (8) follows the identical procedure given in the 
two previous derivations. That is, contributions to the local flux at a point 
Sg in the plane of the open end are each determined and then integrated to give 
the total contribution. In this case the contributing sources are: the orifice, 
the rear wall, and the side wall. These contributions, when set up in integral 
form, become 


n 2 (s 2 ) 





n 3 (x)R(L - x) 



R - So cos p 


+ (L-x)^ - 2Rs 2 cos pj 


dp dx (B18) 


r 0 

After the P integration, equation (B18) becomes 


n 2 ( s 2 ) = n-]J(s 2 ) + 


f 

/h/2 


R /*L 

n 4 (r)K(r, s 2 )dr + / h 3 (x)M(x, s 2 )dx 

•'o 


( 8 ) 


where 
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2 



K(r, s 2 ) 



M(x, s 2 ) = 2R 2 

[(, 


(L - 

- x) 3 + ( 

'f - S |)(L - 

■ x) 



x) 4 + 21 

( b 2 - 4) 

)(L - x) 2 + 1 

[« 2 - 4) 

21 

1 
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APPENDIX C 


DERIVATION OF ANALYTIC EXPRESSION FOR SPECULAR-REFLECTION ANALYSIS 

In equation (12), an expression for the arrival rate passing through a dif- 
ferential hand on a hemisphere is given as 

d(n^) = rij_2 cos 0 ^ sin 9 ^ d0^ (Cl) 

Consider an annulus located at radius s in a plane transverse to the axis 
of the basic model, as shown in sketch (j) . 



U) 

Similar to the method given in the text, particles leaving the orifice in a 
direction 0 are represented by a ray in the two-dimensional plane. Obviously, 

an infinite number of rays can pass through the annulus shown in sketch (j). The 
flux through the annulus is given by the summation of arrival rates of all these 
rays, divided by the area of the annulus: 


n 2 (S2) 2jrs(As) 

To establish the series given in equation (C2), consider the three possible 
paths that rays may take to arrive at the annulus. These are shown typically in 
sketch (j). 

For the first ray 0^_, 




2 cos 0 i sin 0 i (A0 i ) 
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i iiraiim 


II VI II I III 


iii ii i nii mi 1 1 mi ii mi i ii mi rniiNi 



as-l = 


(As)cos 6 1 ^ (Aa) L 

Vs 2 + L 2 


2 , t 2 
s + L 


The next ray © 2 first strikes the side wall and is reflected through the 
annulus. For this ray, 


cos = 


V(2R - s ) 2 + L 2 

(As) cos © 2 


sin ©2 - 


2R - s 


A0 2 = 


V( 2 R - s ) 2 + L 2 
(Ab)L 


V(2E - s > 2 + I? (2S - S)2 + 12 


The third ray, © 3 , first strikes the side wall and is reflected to the axis. 
It then re-reflects from the axis and passes through the annulus. For this ray, 


cos ©j = 


V(2R + s ) 2 + L 2 
(As)cos ©3 


sin 0 j = 


2R + s 


A0 3 = 


V( 2 R + s ) 2 + L 2 
(As)L 


V(2R + s ) 2 + L 2 < 2E + 5)2 + L2 

Equation (C2) can thus he written as 


n ? (s) 


s(As)L 2 (2R - s ) ( As ) L 2 + (2R + s)(As)L 2 


n. 


rts(As) |^2 + l 2 ) 2 j^( 2R - s ) 2 + L 2 ] ^ [(2R + s ) 2 + L 2 ]' 

Equation (C3) reduces to 


+ . 


(C3) 




r 

00 


00 


n 2 (s) 

= ± < 

y 

2oR + s 

y 

2aR - s 

»1 

rts 

L 

2 

[(2oR + s ) 2 + L 2 ] 


[(2aR - s ) 2 + L 2 ] 




k. a= 


(C4) 


Normalized with respect to R and with n^ = h^(jr/4) ( d^/2) , equation (C4) 
may he written to express the relative flux ratio: 
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(C5) 

Inspection of equation (C 5) shows that the term "by term contributions are 
of the order of l/a 3 , so that only a few terms are needed for accurate solution. 
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TABLE I. - VALUES OF RATIO TEST 


AND NUMBER OF ITERATIONS 
OF EQUATIONS (6) AND (7) 


Length to 
radius 
ratio, 

l/r 

Ratio® 

Number of 
iterations 

1 

0.992 • 

20 

2 

. 987 

20 

3 

.981 

20 

4 

. 977 

80 

5 

. 971 

60 

10 

. 968 

20 


a The ratio is defined as fol- 
lows: 


Ratio 



[n 2 (s)/n 1 ]s ds 


1/2 ( 1 ^ 2) 2 


TABLE II. - END-POINT VALUES OF RELATIVE FLUX DISTRIBUTION 
CURVES OF FIGURES 3 to 5 FOR DIFFUSE-REFLECTION 
ANALYSIS OF MODEL SHOWN IN FIGURE 2 


Length to 

Rear wall, 

Side 

wall, 

radius 

n 4 (r 

)/»l 

n 3 (x)/ni 

ratio, 

(fig 

3) 

(fig. 5) 

l/r 









Position 


Radial, r/R 

Axial, x/R 


d2_/2R 

l a 

o a 

(l/r) b 

1 

0.069X10“ 2 

0.029x10-2 

0.043x10-2 

0. 114x10-2 

2 

.178 

.076 

.098 

.108 

3 

.276 

.123 

.148 

.107 

4 

. 371 

.168 

.197 

.102 

5 

. 461 

.213 

. 242 

.101 

10 

. 902 

.429 

. 462 

.098 


Corresponding points. 
"^Corresponding points. 


Open end, 
n 2 (s)/ni 
(fig. 4) 


Radial, s /R 


0 


0. 108X10- 2 
. 131 
.135 
.136 
.136 
.138 


0. 333x10' 
.254 
. 249 
. 248 
. 246 
. 244 
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Figure 1. - Electron-bombardment ion thrustor. 



1 Orifice 

2 Open end 

3 Side wall 

4 Rear wall 


Figure 2. - Basic model for diffuse- and specular-reflection analyses 












Relative flux, r^s)/^- 


Length to radius ratio 



Infinite tube distribu 
tion (ref. 1) normal- 
ized to value of 
L/R = 10 at axis 


Normalized radial distance, s/R 


Figure 4. - Relative flux distribution along 
open end of model shown in figure 2 for 
various length to radius ratios. Diffuse- 
reflection analysis. Radius of orifice, 

0. 0476 units. 







Normalized radial distance, s/R 

(a) Various length to radius ratios. 
Compare with figure 4. 

Figure 6. - Relative flux distribution along 
open end of model shown in figure 2. 
Specular-reflection analysis. 




Relative flux, ng( s) /n 







1 

— - 1 


Approximate solu- 
tion eq„ (13) \ 


-Analytic solu- 
tion eq. (C5) 

- ! I —I - 


Normalized radial distance, s/R 

(t>) Comparison of approximate and analytic 
results for a length to radius ratio 
of 5. 


Figure 6. - Concluded. Relative flux distribu 
tion along open end of model shown in fig- 
ure 2. Specular -re flection analysis. 





Normalized axial distance, x/R 


Figure 7. - Comparison of relative flux distribu- 
tion along side wall of model shown in figure 2 
Length to radius ratio, 4. Diffuse- and spec- 
ular-reflection analysis. 






Model C 


Figure 8. - Section views of ion-thrustor models 
and mesh overlay for specular reflection anal- 
ysis. Length, 4; radius, 1. Normalized units 
Nine equally spaced concentric screen openings 
0.08 units each. 







Mesh overlay: x-inerement 0.4 units; 

r-increment 0.04 units 

Figure 8. - Concluded. Section views of ion-thrustor 
models and mesh overlay for specular reflection 
analysis. Length, 4; radius, 1. Normalized units. 
Nine equally spaced concentric screen openings, 0.08 
units each. 
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Normalized radial distance, r/R 
(h) Normalized axial distance from orifice, 2.8 units. 

Figure 9. - Continued. Specular reflection relative radial 
flux distributions at various normalized axial distances 
from orifice for models shown in figure 8. 
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Normalized radial distance, r/R 

(c) Normalized axial distance from orifice, 3.2 units. 

Figure 9. - Continued. Specular reflection relative radial 
flux distributions at various normalized axial distances 
from orifice for models shown in figure 8. 









Normalized radial di stance , r/R 
(d) Normalized axial distance from orifice, 3.6 units. 

Figure 9. - Concluded. Specular reflection relative radial 
flux distributions at various normalized distances from 
orifice for models shown in figure 8. 








